On maturation of crack patterns 



cn 
o 
o 

> 
O 

o 

(N 



E. A. Jagla 

The Abdus Salam International Centre for Theoretical Physics 
Strada Costiera 11, (34OI4) Trieste, Italy 

Superficial (two dimensional) crack patterns appear when a thin layer of material elastically attached 
to a substrate contracts. We study numerically the maturation process undergone by these crack 
patterns when they are allowed to adapt in order to reduce its energy. The process models the 
evolution in depth of cracks in geological formations and in starch samples ('columnar jointing'), 
and also the time evolution (over thousands of years) of crack patterns in frozen soils. We observe 
an evolution towards a polygonal pattern that consist of a fixed distribution of polygons with mainly 
five, six and seven sides. They compare very well with known experimental examples. The evolution 
of one of these 'mature' patterns upon reduction of the degree of contraction is also considered. We 
find that the pattern adapts by closing some of the cracks and rearranging those in the immediate 
neighborhood. This produces a change of the mean size of the polygons, but remarkably no changes 
of the statistical properties of the pattern. Comparison with the same behavior recently observed 
in starch samples is presented. 
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I. INTRODUCTION 

Consider a thin layer of a solid material elastically at- 
tached to a substrate. If the material contracts (or the 
substrate expands), elastic stresses appear in it. When 
these stresses are sufficiently high, cracks can appear in 
the material, giving rise to a fragmentation process. Well 
known examples of this phenomenon are cracking on mud 
and paints. In these cases the water evaporation pro- 
duces the contraction of the material that is responsible 
for cracking. In other cases, as in the cracking of ce- 
ramic coatings, it is typically the contraction upon cool- 
ing that generates the same phenomenon. Fragmentation 
is known to produce a two dimensional pattern of cracks 
whose statistical properties have been studied theoreti- 
cally [1], experimentally [2] and numerically [3,4]. With 
some variations depending on the particular case, these 
crack patterns are hierarchical structures, with younger 
cracks meeting older ones perpendicularly. Then most 
crack joints are 'T' shaped [5], with the horizontal part 
being older than the vertical part. 

There is however a small number of remarkable cases 
in which fragmentation crack patterns undergo a 'mat- 
uration' process. This means that starting from a hi- 
erarchical pattern as described above, cracks can adapt 
smoothly to optimize its configuration. This optimiza- 
tion process is driven by the tendency of the crack pat- 
tern to reduce its mechanical (elastic plus crack) energy. 
Special conditions have to be fulfilled for this maturation 
to take place. To modify a given crack pattern, cracks 
should be able to displace laterally, and this implies typ- 
ically the surmounting of enormous energy barriers (al- 
though the final state has lower energy than the original 
one). Particular conditions make this lateral displace- 
ment possible in (at least) two different cases. 

One is the case of crack patterns formed on the ground 



of very cold regions of the earth [6], and also in other 
planets [7]. In this case the frozen ground (named 'per- 
mafrost') cracks when the rapidly fallen temperatures of 
winter make the surface contract with respect to lower 
parts of the terrain. This first crack pattern is of the 
kind described above. The cracks get filled with new ice 
and debris, and when temperature rises after winter the 
cracks tend to close. However, the new material that 
filled the cracks is weaker than the old permafrost, and 
the next year cracks open almost on top of the 'scars' of 
first year cracks. However, small lateral variations can 
occur from one year to the next. There are many reasons 
that can make a crack to be shifted laterally in one di- 
rection or the other, from one year to the next. Most of 
these reasons (as for instance inhomogeneities in the ma- 
terials) are not expected to bias the shift of the crack in 
one particular direction. But there is at least one reason 
for a crack to shift in a particular direction, and that is 
the tendency to reduce the energy of the crack pattern. 
In fact, from a statistical point of view it is reasonable to 
expect the crack pattern to adapt in order to reduce its 
energy. This tendency provides a bias for the evolution 
of the crack pattern in permafrost that over thousands of 
years is able to qualitatively modify its appearance [6]. 
In fact, after maturation, crack joints become more 'Y' 
shaped, as this form has lower energy than the 'T' shaped 
original joints. 

The second, better known and more remarkable exam- 
ple of crack pattern maturation takes place in the case of 
columnar jointing. It occurs in basaltic rocks when they 
cool after its expulsion in a volcanic event [8], and also 
in desiccating starch [9] driven by the shrinkage due to 
humidity loss. In both realizations, a superficial pattern 
of cracks very much like the one described in the first 
paragraph first develops in the material. But here, this 
crack pattern penetrates the material as deeper parts of it 
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cool (or dcssicatc). It is this progression into the interior 
that allows the maturation of the pattern to take place, 
now as a function of depth, reaching a polygonal struc- 
ture whose further advance defines prismatic columns. In 
this case there is no true lateral movement of the cracks. 
But a description in terms of lateral movement can be 
given if we choose a reference system that moves with 
the penetrating crack front. 

We use here a recently developed model of fracture [10] 
to describe crack patterns in a two dimensional material 
elastically coupled to a substrate. In the original formu- 
lation of this model cracks have to be pinned in some 
way in order to avoid them to move laterally (since typi- 
cally this movement is unphysical). Here instead, we take 
advantage of this movement (driven by the tendency to 
minimize the energy of the system) to observe how an 
originally disordered pattern becomes polygonal during 
its maturation. We also investigate the way in which a 
stable polygonal pattern is modified when the degree of 
contraction is modified. We observe that some individ- 
ual cracks disappear (terminate, in the 3D language of 
columnar jointing) when contraction is reduced, giving 
rise to local rearrangements in the pattern. This mech- 
anism provides a way to change the mean width of the 
columns as a function of depth in the basalt formations 
and in starch, and it has been observed to occur in this 
last case. We finish with a discussion on what the typical 
width of columns in three dimensional formations is. 



II. THE NUMERICAL MODEL 

We use a technique recently developed [10] to treat 
fracture and cracks in the context of phase field mod- 
eling [11]. The free energy of the system is written in 



terms of the strain tensor e.; 



\/2{dui/dxj + duj /dxi), 



with u(r) being the local displacement field. We choose 
the form of the free energy in such a way that it re- 
duces to the normal elastic energy for small strains, but 
for large strains it is able to describe cracks. This is 
achieved by a saturation of the free energy of the system 
for large values of Sij . The inclusion in the free energy of 
terms proportional to gradients of e produces a smooth- 
ing of cracks, which although artificial, is however very 
important to us. On one hand it makes the description 
isoptropic and insensitive to the numerical mesh we use 
in the calculation (as long as the discretization is much 
thinner than the smoothing distance of the fracture) . On 
the other hand it allows the cracks (that in the regular- 
ized theory could be pictorially described as 'solitons') to 
move around the system to find configurations of lower 
energy. This wandering will model the maturation of the 
crack pattern. The free energy is taken to be rotation- 
ally invariant, in order to describe cracks in an isotropic 
material. 



To write down explicitly the equations we actually 
solved in our two dimensional geometry, we first intro- 
duce the following notation for the independent compo- 
nents of e [12,13] 



ei = (£11 +£22)/2 

62 = (£11 - £22)/2 

63 = £12 = £21 



(1) 



which are named respectively the dilation, dcviatoric and 
shear components. These three variables are not inde- 
pendent. They satisfy the St. Venant compatibility con- 
straint [12,13] 

[dl + 5^)ei - [dl - dl)e2 - 2d,dyes = 0. (2) 

The free energy density is 



F{e) = 



[l + FO{e)/fO] 



(3) 



where 



F°{e) = B{e, e^f + ^^ ((e^ - e^f + {e, - el)') (4) 

and B and /it are related to the two dimensional bulk 
and shear modulus of the material. e^(r) are externally 
controlled functions that allow to prescribe the locally 
preferred state of the system, and g{r) is another (pos- 
itive) function that will be used to model some random 
inhomogencities in the system. The limiting value /o of 

for £ — > 00 (assuming g = 1), is related to the crack 
energy in the model. 

Regularization of cracks is provided by a gradient term 
Fg in the free energy density, that we choose to be of the 
form 



Fg= J2 ^ii'^^if 



(5) 



=1,2,3 



where we have to choose 0:2 = as to retain rotational 

invariance. 

An additional ingredient that has to be added here 
with respect to the basic model of Ref. [10] is the inclusion 
of the clastic energy density F^i of the system attached 
to the substrate. In terms of the displacement variables 
u, this elastic energy can be written in the form 



dr'Fei=j / dr2|u(r)| 



(6) 



where 7 measures the stiffness of the interaction with 
the substrate. As we take the components of s to be our 
basic variables, we have to recast this energy in terms of 
them. This can be easily done in the Fourier space, and 
the result is 



|g2(k)|^ + |g3(k)|^ 
fc2 



(7) 
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Where ei(k) are the Fourier transforms of the original 
ej(r). The equations of motion are taken to be of the 
over damped form, namely 



where 



dei{v) SF 



F = Jdv'' (F + Fei+Fg) 



(8) 



(9) 



The Saint Venant constraint is implemented by using a 
Lagrange multiplier. 



III. RESULTS 

We did the simulations on a square mesh of 512x512 
elements, using periodic boundary conditions. Starting 
from the flat configuration ei(r) = 62(1") = e3(r) = we 
simulated a uniform and abrupt contraction of the system 
by taking e?(r) = c, e^(r) = 0, ei(r) = 0, c = 0.7. We 
introduce also a finite disorder, taking g(r) in (3) to be a 
random function on the lattice, uniformly distributed be- 
tween 0.75 and 1.25. We keep B and the mesh discretiza- 
tion 5x as scale-fixing parameters, and take /i = 0.5i3, 
ai = 0.5Bdx^ {i = 1, 2, 3), 7 = 0.0025B/(5a;2, and 
/o = 0.5-8. The crack energy per unit length rj is then 
77 = faSx. Under these conditions we solve the evolution 
equations. We see in Figs. 1 and 2 snapshots of the time 
evolution of the system. The figures are done by marking 
the points in which ei > 0.5. According to our definition 
of the free energy (3) and (4), this is a reasonable defini- 
tion of 'broken' elements (results do not depends strongly 
on the threshold value used) . In this way we arc basically 
plotting the cracks present in the system. It is important 
to note that due to the finite value of ai , broken elements 
do not form strictly onc-dimcnsional 'strings' in the sys- 
tem, but they clusterize, making cracks acquire a finite 
width, as it is apparent in the pictures. This is a crucial 
point to simulate an isotropic system. 

We can distinguish two different stages in the tem- 
poral evolution. During the nucleation stage (Fig. 1) 
cracks appear rather disorderly in the system and prop- 
agate around. The pattern that forms is very dependent 
on many details of the simulation, as for instance the 
amount of disorder present. This is the kind of pattern 
we have described in the introduction as a fragmenta- 
tion pattern. During a second stage the maturation of 
the pattern occurs (Fig. 2). This is observed as a pro- 
gressive lateral displacement of the cracks towards a con- 
figuration of lower energy. It is necessary to emphasize 
again that in standard fragmentation processes this mat- 
uration cannot take place, as cracks are rigidly located in 
their positions. In our numerical model cracks can in fact 
move laterally, since this does not imply the surmounting 
of a large energy barrier. 



The lateral movement of cracks in our model is however 
rather slow compared to its nucleation, and that is why it 
is not seen on the timescale of the nucleation stage. We 
stress that we are not forcing the crack pattern to become 
polygonal, or cracks to terminate onto other cracks, it is 
the system itself that prefers this kind of configuration as 
this reduces its energy. The final, stable pattern is that 
at the bottom right of Fig. 2. It corresponds to a rela- 
tive minimum of the energy of the system, the absolute 
minimum being a perfect hexagonal pattern with a poly- 
gon size (calculated numerically with the same model) as 
indicated also on Fig. 2. The mature pattern contains 
mostly polygons of five, six, and seven sides, and a small 
number with four, and eight sides. They are statistically 
very similar to those in real columnar formations (see Fig. 
4 below and Fig. 8 in Ref. [14]). We note that the mean 
area of polygons for different number of sides follows a 
linear relation, known as the Lewis law, after he encoun- 
tered it in other two dimensional patterns [15]. This law 
follows if the pattern is assumed to be maximally random 
[16]. 

The present results can be compared with those ob- 
tained previously [14] using a phenomenological model 
for the energy of the cracked material. The present ap- 
proach is however much more general than that in Ref. 
[14]. Here, we are not assuming any phenomenological 
form of the energy as a function of the areas of the poly- 
gons, the energy of the system builds up from the free 
energy presented in the previous section. In addition, 
crack segments are not forced here to be straight, and in 
fact we can see in the last panel of Fig. 2 that some of 
them are slightly curved. The curvature occurs particu- 
larly when there is a large difference between the areas of 
polygons on both sides of the crack segment, always curv- 
ing it in the direction in which areas tend to be closer. 
The reason for this is again energetic: slightly curving 
a crack does not pay much crack energy, but produces 
a gain in elastic energy if the areas of the two adjacent 
polygons tend to become closer to each other. This cur- 
vature has been in fact observed to occur in a full three 
dimensional calculation for a simple geometry [17]. 

An interesting problem to be investigated with the 
present model is the way in which a stable polygonal pat- 
tern changes when there are changes in the parameters 
that control the extent of contraction. As an outcome 
of this analysis we will get an idea of the expected evo- 
lution of the patterns down in the columnar formation 
(after the first maturation), since the thermal stresses in 
deeper parts of the material are lower than close to the 
surface. Since a lower grade of contraction corresponds 
to an ideal pattern with larger polygons, we may wonder 
what is the way (if any) in which one of our patterns 
adapts to the new conditions. We present in Fig. 3 the 
results of simulations when the extent on contraction c 
is reduced. We see that there is an increase in the mean 
area of the polygons when c is reduced. The area in- 
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crease is not homogeneous over all polygons, but occurs 
due to the disappearance of particular crack segments, 
merging two (or three) adjacent polygons into one. After 
the disappearance of the crack segments there is a local 
rearrangement of the pattern which adjusts to the new 
configurations. Those regions in which no crack disap- 
pear remain perfectly stable despite the change in the 
contraction. 

The evolution of the mean area A of the polygons as 

a function of c is shown in Fig. 4 along with the ideal 
area A*'' of the hexagons in the perfect hexagonal pattern 
of minimum energy. We see that the evolution tends to 
follow that of the ideal structure, although A is always 
smaller than A^'^ [18]. We note also that in the present 
model there is a critical value of c 0.42) for which ^4"^ 
diverges, and wc expect the same occurs for non ideal 
patterns. This happens because the elastic energy per 
unit area gained when generating a polygonal crack pat- 
tern decays very rapidly when the size of the hexagons 
increases sufficiently. The sum of this elastic energy plus 
fracture energy may not have a minimum with respect to 
the area of polygons if the degree of contraction c is too 
small. Note that the same does not apply to a real three 
dimensional columnar case (see next section). 

It is remarkable that the statistical properties of the 
pattern do not change appreciable during this relaxation 
stage. In Fig. 5 we see that despite a change in the 
mean area by more than fifty percent, statistical distri- 
bution of polygons by number of sides and areas remain 
constant within numerical fluctuations associated to the 
finite size of the system [19]. Note that it is precisely 
the 'imperfection' of the crack pattern that makes pos- 
sible the adaptation of the mean area to a condition of 
lower contraction. For a perfect hexagonal pattern it 
is impossible to find a way to adapt the pattern slightly 
and obtain another hexagonal pattern with slightly larger 
polygon area. In our case the mean area of the pattern is 
increased by making some crack segments between poly- 
gons disappear. 

IV. THE PROBLEM OF THE TYPICAL COLUMN 
WIDTH IN COLUMNAR FORMATIONS 

The two dimensional model we have studied is per- 
fectly well defined, and provides values for the size of the 
polygons in the ideal hexagonal pattern that minimizes 
the total energy of the system. We want to comment at 
this point to what extent these two dimensional results 
can be applied to the full three dimensional columnar 
problem. For a straightforward application to be pos- 
sible, the elastic energy stored in a columnarly cracked 
three dimensional material should be stored in a layer 
around the crack front of thickness w, in such a way 
that w is much smaller than the typical column width I 
(Z ~ A^/^). If this is satisfied, the two dimensional de- 



scription is directly applicable. The only consideration 
to be made is that two-dimensional variables have to be 
scaled from the three dimensional variables using w. For 
instance, the effective crack energy per unit length rj and 
elastic constants B and fi of the two dimensional descrip- 
tion are obtained from the real three dimensional values 
as wri^^'^\ wB'^^'^^ and w/i'-^'^h Unfortunately the condi- 
tion w « I is never satisfied. In fact, the elastic energy 
of a columnar formation is stored in a portion of thick- 
ness w ^ I around the crack front [20] . The coincidence of 
the statistical properties of our two dimensional patterns 
and those in true three dimensional cases indicates that 
these properties are robust with respect to this difference. 
However, the calculation of the ideal size of the perfect 
hexagonal pattern (and then an estimation of the typical 
size of non-perfect real patterns) has to be reconsidered 
for the three dimensional case. In fact, in our two dimen- 
sional model, in which a layer of material is attached to 
a substrate, there is an ideal hexagonal pattern of well 
defined polygon size that minimizes the energy of the 
system. The application of the same principle of min- 
imizing the total energy leads in the three dimensional 
case to nonsense: the contribution of the fracture energy 
to the total energy is always much larger than the elastic 
contribution. A minimum can only be obtained with no 
cracks at all. 

The correct way to pose the problem of the typical size 
of the polygons in three dimensions is the following: in 
the three dimensional case, we have a temperature profile 
that we assume to be dependent only on depth z, passing 
more or less steeply from Tq at z ^ — oo to Ti > Tq at 
z — > oo, in such a way that a temperature front can be 
defined. At time t = the temperature front is assumed 
to be located at z = 0. We will consider the idealized 
case in which the temperature profile is rigidly displaced 
towards the interior as a function of time, with some 
fixed velocity w, namely T{z. t) = T{z — vt). We assume 
that a stable polygonal pattern of fractures has formed, 
and that its front is located at some depth zq, which 
moves down locked to the temperature profile, namely 
Zq = Zr + vt, whcrc Zr mcasurcs the relative position of 
the crack pattern and the temperature profile. The value 
of Zr depends mainly on the typical size of the pattern I 
and the overall temperature difference AT = Ti — Tq. A 
previous stability analysis has shown in a simplified case 
[20] that under the present conditions, patterns with dif- 
ferent I can propagate in a stable manner, with Zr being a 
decreasing function of I, namely larger patterns are more 
retarded with respect to the temperature front. However 
there is a limit to this stable propagation. If I or AT 
are too small, the crack front becomes unstable: not all 
cracks can propagate. It is tempting to argue (and this 
is also based upon what is observed in three dimensional 
starch samples, see below) than in this case some crack 
segment will remain halted, and the rest of the pattern 
propagates. In this way / is effectively increased and the 
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crack front becomes closer to the temperature front, in 
such a way that the new pattern is now stable. In a situ- 
ation in which AT decays smoothly with time (whereas 
at the same time the front penetrates the material), we 
may expect that the pattern will always be located at the 
value Zcr that marks the limit between stable and unsta- 
ble propagation. This is the condition that determines 
the size of the columns in terms of the temperature pro- 
file, the elastic properties of the material and the crack 
energy. For the case of a sharp temperature jump and 
generalizing the two dimensional expressions for elastic 
and fracture energy in Ref. [20] , we obtain that the crack 
front is located precisely at the border between stable 
and unstable regions when B{aAT)^l/r] is some constant 
value k of order unity (this value is not easy to calcu- 
late). Here a is the thermal expansion coefficient and B 
is a typical elastic constant of the material. From here 
we obtain the typical width of the columns as 

l = k^{aAT)-^ (10) 

The typical size is then positively correlated to the crack 
energy, and negatively correlated with the elastic stiffness 
of the material, both facts being qualitatively reasonable. 
The size I is also proportional to the negative second 
power of the temperature jump responsible for cracking. 
We should keep in mind however that this result is valid 
only for the assumed sharp step form of the temperature 
profile. In other cases we should search for the critical 
position of the crack front Zcr along the lines used in 
Ref. [20]. Note that a AT plays in three dimensions the 
role of the degree of contraction c in our simulations. 
The preceding formula indicates that the typical size of 
polygons diverges only when aAT 0, contrary to the 
critical value of c we found in two dimensions (see Fig. 4). 
This would indicate that if the driving force for cracking 
is slowly reduced when going deeply into the material (as 
may occur be due to the higher difficulty to expel heat -or 
humidity in starch- through the upper material) the size 
of the pattern should adapt by increasing their typical 
size, but it would never stop abruptly. 

Recent tomography experiments in starch samples [21] 
show that termination and rearrangement of cracks seem 
to be in fact the main mechanism by which the polygonal 
pattern evolves in depth. In starch samples the humidity 
gradients are expected to be reduced when going deeper 
into the sample, and that is why the typical width of the 
columns tends to increase. However, a quantitative veri- 
fication of a relation like (10) (or the equivalent one for a 
more realistic time dependent temperature or humidity 
profile) is not possible at present as it would require the 
in situ determination of the temperature profile under 
which the cracks form, and not only the observation a 
posteriori of column thickness as a function of depth. 



V. CONCLUSIONS 

We have studied numerically the formation and matu- 
ration process of a two dimensional crack pattern that is 
allowed to adapt to find configurations of minimum en- 
ergy. The original cracks appear in a rather disordered 
way, but the pattern naturally evolves towards a polygo- 
nal configuration with well defined statistical properties. 
We argue that this maturation process occurs in cracks 
patterns on the ground of arctic regions (permafrost) 
and effectively in the columnar jointing of basalts and 
starches, as a function of depth. Our model allows also 
to study the evolution of mature polygonal patterns when 
the extent of contraction is reduced. We have found that 
in this case the pattern adapts by closing ('terminating' 
in the three dimensional interpretation) some cracks and 
rearranging those cracks in the immediate neighborhood. 
This evolution has been recently observed to occur in 
starch samples. Although it does not contain all features 
of the full three dimensional problem, our approach pro- 
duces patterns of very good statistical agreement with 
real ones. The issue of the typical scale of the three 
dimensional pattern is beyond the reach of the two di- 
mensional model, and we have provided for this case a 
plausible description that relates the typical size of the 
polygons with the elastic and thermal properties of the 
material, and with the details of the temperature profile. 



[1] K. Leung and J. V. Andersen, Europhys. Lett. 38, 589 

(1997) ; S. Kitsunezaki, Phys. Rev. E, 60, 6449 (1999). 
[2] A. Groisman and E. Kaplan, Europhys. Lett. 25, 415 

(1994); P. Mcakin, Science 252, 226 (1991); W. Korneta, 
S. K. Mondiratta, and J. Menteiro, Phys. Rev. E 57, 3142 

(1998) ; K. A. Shorlin, J. R. de Bruyn, M. Grahan, and 
S. W. Morris, Phys. Rev. E 61, 6950 (2000). 

[3] T. Hornig, I. M. Sokolov, and A. Blumen, Phys. Rev. E 

54, 4293 (1996). 
[4] K. Leung and Z. Neda, Phys. Rev. Lett. 85, 662 (2000). 
[5] Except those occurring as a consequence of bifurcation 

of single cracks. 
[6] A. H. Lachenbruch, Geol. Soc. Am. Spec. Pap. 70, 69 

(1962); R. S. Slette, B. Hallct, and R. C. Fletcher, J. 

Gcophys. Res. 108(E4), 8044 (2003). 
[7] H. Hiesinger and J. W. Head, J. Geophys. Res. 105(E5), 

11999 (2000). 

[8] J. Walker, Sci. Am. 255, 178 (1986); D. L. Peck and T. 
Minakami, Geol. Soc. Am. Buh. 79, 1151 (1968); D. Wea- 
rie and C. O'CarroU, Nature (London) 302, 240 (1983); 
P. E. Long and B. J. Wood, Geol. Soc. Am. Bull. 97, 
1144 (1986); J. M. DeGraff and A. Aydin, Geol. Soc. 
Am. Bull. 99, 605 (1987); A. Aydin and J. M. Dcgraff, 
Science 239, 471 (1988); P. Budkewitsch, P. Robin, J. 
Volcanol. Geotherm. Res. 59, 219 (1994). 



5 



[9] G. Miillcr, J. Geophys. Res. 103, 15239 (1998); G. 

Miiller, J. Struct. Geol. 23, 45 (2001) 
[10] V. I. Marconi and E. A. Jagla, to bo published. 
[11] H. Emmerich, The Diffuse Interface Approach in Mate- 
rials Science, (Springer, Berlin, 2003). 
[12] S. Kartha, J. A. KrumhansI, J. P. Sethna, and L. K. 

Wickham, Phys. Rev. B 52, 803 (1995). 
[13] S. R. Shenoy, T. Lookman, A. Saxena, and A. R. Bishop, 

Phys. Rev. B 60, R12537 (1999); T. Lookman, S. R. 

Shenoy, K. O. Rasmussen, A. Saxena, and A. R. Bishop, 

Phys. Rev. B 67, 024114 (2003). 
[14] E. A. Jagla and A. G. Rojo, Phys. Rev. E 65, 026203 

(2002). 

[15] F. T. Lewis, Anat. Rec. 38, 341 (1928); ihid, Am. J. Bot. 
31, 619 (1944). 

[16] D. Weaire and N Rivier, Contemp. Phys. 25, 59 (1984); 

N. Rivier and A. Lissowski, J. Phys. A 15, L143 (1982). 
[17] R. Saliba and E. A. Jagla, J. Geophys. Res. 108(B10), 

2475 (2003). 

[18] This difference can be reasonably considered to be a con- 
sequence of the metastabilities in the process we are sim- 
ulating. In fact, if at some point we revert the evolution 
of contraction, making now c increase, we have observed 
the nucleation of new crack segments, typically bisecting 
the largest polygons. In this evolution the values of A 
we get are larger then ^4*'', namely the evolution of A as 
a function of c has a strong hysteresis that encloses the 
ideal value corresponding to the ideal energy minimum. 

[19] In Fig. 5(a), it seems that the first stable pattern (with 
the original value of contraction c = 0.7) is slightly differ- 
ent to the rest, having less polygons with six, and more 
with five and 7 sides. We cannot tell at this point if this 
is systematic or associated to statistical fluctuations. 

[20] E. A. Jagla, Phys. Rev. E 65, 046147 (2002). 

[21] L. Goehring and S. W. Morris, Bulletin of the American 
Physical Society, 48, 611 (2003); ibid, private communi- 
cation. 




t=1.5 T i-Z T 



FIG. 1. Appearance of a typical fragmentation pattern dur- 
ing the first stage of the evolution. The time scale t is given 
by = XB. The contraction imposed is c = 0.7 




FIG. 2. Maturation of the fragmentation pattern at longer 
times (note the change in of time intervals with respect to 
previous figure). The lateral displacements of cracks allows 
the system to reach a stable state (bottom right) which is a 
local energy minimum. The size of the hexagon in the per- 
fect pattern that corresponds to the absolute minimum of the 
energy is indicated. 
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FIG. 3. Evolution of the mature pattern of Fig. 2 (bottom 
right) upon reduction of the extent of contraction c. Note the 
disappearance of some cracks and the local rearrangement 
that occur. To facilitate the visualization we plot also the 
immediately previous pattern. For the first panel the previous 
pattern is the last one in Fig. 2. 
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FIG. 5. (a) Frequency of polygons with different number 
of sides and (b) mean area of polygons with different number 
of sides (normalized to the mean area of all polygons) for 
patterns obtained by reducing c. 
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FIG. 4. Evolution of the mean area A of the polygons as 
a function of the degree of contraction c, when c is reduced 
from larger to smaller values, and ideal value A!' of the area of 
polygons in the perfect hexagonal pattern that minimizes the 
energy of the system. There is a critical value of contraction 
(c ~ 0.42) below which the uncracked configuration is the one 
with minimum energy. 
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